# Purpose: Replication file for the paper "Networks, Dyads, and the Social Relations Model" 
# Author: Cassy Dorff & Michael D. Ward

#####################################################
# set up
#####################################################
rm(list=ls())
datapath<-c("/Users/cassydorff/ProjectsGit/boulder/code/Replication Files/")
setwd(datapath)
library(foreign)
library(network)
library(qgraph)
library(foreign)
library(MASS)
library(xtable)

#####################################################
# load all the data
#####################################################

# TRADE DATA for section 3 illustration
load("data/trade.rda")


# MID DATA
load("data/Soc.Matrix.dta")  #a list of MIDs matrices for section 4.1
load("data/Soc.Matrix.Agg.acc.dta")   #aggregated matrix of counts MIDs for 4.1 
load("data/Soc.Matrix.Agg.dta")   #matrix of binary onsets MIDs for 4.1 

# TRADE DATA (trade matrix for GBME and SRM section 4.2)
dd<-read.csv(paste(datapath,"data/exIMF2000.csv",sep="/"), header=T)

# MORROW DATA (and network creation for section 4.3)
### data and matrices for last time slice
morrow.pows.t5.dat<-read.dta("data/morrow_pows_t5.dta")
morrow.pows.t5.dem<-read.dta("data/morrow_pows_t5_dem.dta")
morrow.pow5.edges<-morrow.pows.t5.dat[,1:2]
morrow.pows.t5.ally<-read.dta("data/morrow_pows_t5_ally.dta")
morrow.pow5.ally<-network(morrow.pows.t5.ally[,1:2], directed=F)

set.edge.attribute(morrow.pow5.ally, "ally", morrow.pows.t5.ally$atopally)
morrow.pow5.net<-network(morrow.pow5.edges, directed=T)
set.edge.attribute(morrow.pow5.net, "nocomply", morrow.pows.t5.dat$violator_4_ordinal_comply)
summary.network(morrow.pow5.net)

morrow.pows.t5.ally<-read.dta("data/morrow_pows_t5_ally.dta")
morrow.pow5.matrix<-as.sociomatrix(morrow.pow5.net, attrname="nocomply")
morrow.war5.matrix<-as.sociomatrix(morrow.pow5.net)
morrow.ally5.matrix<-as.sociomatrix(morrow.pow5.ally, attrname="ally")


#####################################################
#Functions
#####################################################

source("functions/SRM.R") 
source("functions/gbme.asym.r")

#####################################################
#Analysis & Graphics
#####################################################

###########################
### TRADE DATA example section 3
###########################
# script for Trade Example in Dyad Paper
# Manufacturing trade among  US, Japan, China, S.Korea, and N. Korea in 2011
# US is 111, JAPAN is 158, CHINA is 924, S. Korea is 542, N. Korea is 954

names<-c("United States","Japan","China","S. Korea","N. Korea")
get.these<-c("111","158","924","542","954")
tr.11<-trade[[2]]
tr.11<-round(log(tr.11+1),2)
tr.11.5<-tr.11[get.these,get.these]
dimnames(tr.11.5)[[1]]<-dimnames(tr.11.5)[[2]]<-names
table1<-tr.11.5

footnote7<-dyads(tr.11.5)

footnote<-cbind(round(footnote7$actor.effect.i,2),round(footnote7$partner.effect.i,2))
library(xtable)
xtable(footnote)


###########################
### MID DATA (section 4.1)
###########################
#  SRM ON binary MATRIX & plots
# run the SRM
out<-dyads(Soc.Matrix.Agg)

ncty<- dim(Soc.Matrix.Agg)[[1]]
sort.ai.agg.bin<-out$actor.effect.i[order(out$actor.effect.i[,1]),]

#figure 2a 
plot(c(0,ncty+5) , c(-.02, .25), type="n", ylab="Actor Effect" , axes=FALSE, ann=FALSE)
axis(1, at=c(1, ncty),labels=c("low", "high" ),1,cex=.13)
title(main="Sender Effects, From Binary Matrix for 50 Years \n Top 10 Senders Indicated")
points(sort.ai.agg.bin, pch=19, cex=.357)
abline(h=0) 
top10<-names(sort.ai.agg.bin[188:197])

lines(c(48.5, ncty-2), c(.206, .195), col="gray90")
text(20, .206, "Yugoslavia/Serbia", cex=.9)

lines(c(28, ncty-3), c(.19, .193), col="gray90")
text(20, .19, "Iraq", cex=.9)

lines(c(32.5, ncty-4), c(.18, .1832), col="gray90")
text(20, .18, "Russia", cex=.9)

lines(c(28, ncty-5), c(.167, .1624), col="gray90")
text(20, .167, "Iran", cex=.9)

lines(c(32.5, ncty-6), c(.148, .137), col="gray90")
text(20, .149, "Jordan", cex=.9)

lines(c(28, ncty-7), c(.1317, .1320), col="gray90")
text(20 , .1320, "USA", cex=.9)

lines(c(32.5, ncty-8), c(.097, .0911), col="gray90")
text(20 , .097, "France", cex=.9)

lines(c(31, ncty-9), c(.083, .08020), col="gray90")
text(20, .083, "China",cex=.9) 

lines(c(27, ncty-10.5), c(.07, .067), col="gray90")
text(20, .07, "UK", cex=.9)

lines(c(30, ncty-10.5), c(.058, .065), col="gray90")
text(20, .058, "Israel",cex=.9)

# Now Target (Figure 2b)
sort.ai.agg.bin<-out$partner.effect.i[order(out$partner.effect.i[,1]),]

plot(c(0,ncty+5) , c(-.02, .25), type="n", ylab="Actor Effect" , axes=FALSE, ann=FALSE)
axis(1, at=c(1, ncty),labels=c("low", "high" ),1,cex=.13)
title(main="Target Effects, From Binary Matrix for 50 Years \n Top 5 Targets Indicated")
points(sort.ai.agg.bin, pch=19, cex=.357)
abline(h=0) 
top10<-names(sort.ai.agg.bin[188:197])

lines(c(50, ncty-2.5), c(.1944, .1944), col="gray90")
text(21, .1944, "Yugoslavia/Serbia", cex=.9)

lines(c(32, ncty-3), c(.135, .128), col="gray90")
text(20, .135, "Russia", cex=.9)

lines(c(29, ncty-4), c(.1221, .1221), col="gray90")
text(20 , .1221, "USA", cex=.9)

lines(c(26, ncty-5), c(.0962, .0962), col="gray90")
text(20, .0962, "UK", cex=.9)

lines(c(28, ncty-5.7), c(.0825, .0825), col="gray90")
text(20, .0825, "Iraq", cex=.9)


# now run the SRM on accumulated counts matrix
out<-dyads(Soc.Matrix.agg.acc)
# Figure 2 (a and b)
ncty<- dim(Soc.Matrix.Agg)[[1]]
sort.ai.agg.on<-out$actor.effect.i[order(out$actor.effect.i[,1]),]
# FIgure 3a 
plot(c(0,ncty+5) , c(-.025, .8),  ylim=c(-.025, .8), type="n", ylab="Actor Effect" , axes=FALSE, ann=FALSE)
axis(1, at=c(1, ncty),labels=c("low", "high" ),1,cex=.13)
title(main="Sender Effects, Aggregated over 50 Years \n Top 10 Senders Indicated")
points(sort.ai.agg.on, pch=19, cex=.357)
abline(h=0) 
top10<-names(sort.ai.agg.on[188:197])

lines(c(32, ncty-1.8), c(.77, .77), col="gray90")
text(20, .776, "Russia", cex=.9)

lines(c(30, ncty-3.5), c(.54, .512), col="gray90")
text(22 , .54, "USA", cex=.9)

lines(c(29, ncty-4), c(.49, .493), col="gray90")
text(22, .49, "Iraq", cex=.9)

lines(c(32, ncty-5), c(.44, .46), col="gray90")
text(22, .44, "China", cex=.9) 

lines(c(30, ncty-6), c(.4, .453), col="gray90")
text(23, .4, "Iran", cex=.9)

lines(c(32, ncty-7), c(.31, .27), col="gray90")
text(23, .31, "Syria", cex=.9)

lines(c(47, ncty-7.5), c(.27, .246), col="gray90")
text(19, .27, "Yugoslavia/Serbia", cex=.9)

lines(c(36, ncty-8.5), c(.22, .235), col="gray90")
text(25, .22, "Turkey", cex=.9)

lines(c(35, ncty-9.8), c(.17, .194), col="gray90")
text(25,.17, "Egypt", cex=.9)

lines(c(34, ncty-11), c(.12, .189), col="gray90")
text(25, .12, "Israel",cex=.9)


# Now Target effects (figure 3b)
sort.ai.agg<-out$partner.effect.i[order(out$partner.effect.i[,1]),]
pdf(file = paste(gd,"50yr_aggT.pdf",sep="/"),width=6, height=5)
 
plot(c(0,ncty+5) , c(-.025, .8) ,ylim=c(-.025,.8), type="n", ylab="Actor Effect" , axes=FALSE, ann=FALSE)
axis(1, at=c(1, ncty),labels=c("low", "high" ),1,cex=.13)
title(main="Target Effects, Aggregated over 50 Years \n Top 10 Targets Indicated")
points(sort.ai.agg, pch=19, cex=.357)
abline(h=0) 
top10<-names(sort.ai.agg[188:197])

lines(c(28, ncty-2), c(.64, .637), col="gray90")
text(20, .637, "USA", cex=.9)

lines(c(32, ncty-3), c(.49, .431), col="gray90")
text(19, .49, "Russia", cex=.9)

lines(c(48, ncty-4), c(.43, .4), col="gray90")
text(20, .43, "Yugoslavia/Serbia", cex=.9)

lines(c(30, ncty-5), c(.38, .365), col="gray90")
text(20, .38, "Israel",cex=.9)

lines(c(27, ncty-6), c(.33, .312), col="gray90")
text(20, .33, "Iraq", cex=.9)

lines(c(30, ncty-7), c(.28, .295), col="gray90")
text(20, .28, "China", cex=.9)

lines(c(25, ncty-8), c(.23, .265), col="gray90")
text(20, .23, "UK", cex=.9)

lines(c(30, ncty-9), c(.18, .255), col="gray90")
text(20, .18, "Japan", cex=.9)

lines(c(28, ncty-10), c(.14, .172), col="gray90")
text(20, .14, "India", cex=.9)

lines(c(30, ncty-11), c(.09, .153), col="gray90")
text(20, .09, "Egypt", cex=.9)

# apply the SRM to each matrix in Soc.Matrix (a list of MID matrices over time)
output.4<-list()
year <- 1816
		for (i in Soc.Matrix){
			print(year)
			output.4<- c(output.4,list(dyads(i)))
			names(output.4)[length(output.4)] <- paste("year",year, sep="_")
			year <- year+1
		}
output.4<-output.4


# pull out all of the individual SRM stats for all the years
actor.effect.i <- lapply(output.4, function(year) year$actor.effect.i)
partner.effect.i <- lapply(output.4, function(year) year$partner.effect.i)
unique.effect.ij <- lapply(output.4, function(year) year$unique.effect.ij) #mat
unique.variance<- lapply(output.4, function(year) year$unique.variance) #1/yr
relationship.covariance<-lapply(output.4, function(year) year$relationship.covariance) #1/yr
actor.variance<- lapply(output.4, function(year) year$actor.variance) #1/yr
partner.variance<-lapply(output.4, function(year) year$partner.variance) #1/yr
actor.partner.covariance<-lapply(output.4, function(year) year$actor.partner.covariance) #1/yr
colmeans <- lapply(output.4, function(year) year$colmeans)


# Figure 4
	country<-c("USA", "GMY", "UKG", "RUS", "IND", "ITA", "FRN", "CHN", "JPN", "AUS")
	ctylist <- data.frame(matrix(NA,length(actor.effect.i),length(country)))
				for (cty in 1:length(country)){
				for ( i in 1:length(actor.effect.i)){
				ctylist[i,cty]<-	actor.effect.i[[i]][match(country[cty], row.names(actor.effect.i[[i]]))]		
				}
				}				
names(ctylist) <- country
years<-1816:2001
row.names(ctylist)<-years
ctylist[ctylist==0]<-NA 
years <- 1950:2001 	

#figure 4a
plot(c(1, length(years)),c(-.05, .3) ,type="n",xlab="Years", ylab="actor effects" , axes=FALSE, ann=FALSE)
axis(1, at=seq(1,length(years),1), lab=c(years))
axis(2, las=1)
title(main="Individual Effect of Actor as Sender")
points(1:length(years), ctylist[135:186,1], col="red", type="l", pch=19, cex=.6,lwd=2) #USA	
points(1:length(years), ctylist[135:186,4], col="slateblue", type="l", pch=19, cex=.45,lwd=2) #RUS 
points(1:length(years), ctylist[135:186,8], col="black", type="l", pch=19, cex=.5,lwd=2) #CHN
legend(40,.18, legend=c("USA","Russia","China"), border=F, bty="n", fill=c("red", "slateblue","black"), inset=.1, ncol=1)


# figure 4b plotting target effects
country<-c("USA", "RUS", "CHN")
ctylist.pe <- data.frame(matrix(NA,length(partner.effect.i),length(country)))
				for (cty in 1:length(country)){
				for ( i in 1:length(partner.effect.i)){
				ctylist.pe[i,cty]<-	partner.effect.i[[i]][match(country[cty], row.names(partner.effect.i[[i]]))]		
				}
				}
				
names(ctylist.pe) <- country
ctylist.pe[ctylist.pe==0]<-NA
plot(c(1, length(years)),c(-.02, .1)  ,type="n",xlab="Years", ylab="actor effects" , axes=FALSE, ann=FALSE)
axis(1, at=seq(1,length(years),1), lab=c(years))
axis(2, las=1)
title(main="Individual Effect of Actor as Target")
points(1:length(years), ctylist.pe[135:186,1], col="red", type="l", pch=19, cex=.6,lwd=2) #USA	
points(1:length(years), ctylist.pe[135:186,2], col="slateblue", type="l", pch=19, cex=.45,lwd=2) #RUS 
points(1:length(years), ctylist.pe[135:186,3], col="black", type="l", pch=19, cex=.5,lwd=2) #CHN
legend(40,.06, legend=c("USA","Russia","China"), border=F, bty="n", fill=c("red", "slateblue","black"), inset=.1, ncol=1)


# figure 5a 
plot(ylim=c(-.1,.25),xlim=c(-.1,.25),bty="n",las=1,ctylist[,1],ctylist[,4],pch=19,xlab="USA Actor Effect",ylab="Russia Actor Effect",axes=T,
     col=c(rep("gray67",123),rep("red",length(124:130)),rep("gray50",length(131:174)),rep("gray90",length(175:186))), cex=.9
)
abline(h=summary(ctylist[,1])[3])
abline(v=summary(ctylist[,4])[3])
arrows(summary(ctylist[,1])[3],summary(ctylist[,4])[3],ctylist[124,1],ctylist[124,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[125,1],ctylist[125,4],ctylist[126,1],ctylist[126,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[124,1],ctylist[124,4],ctylist[125,1],ctylist[125,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[126,1],ctylist[126,4],ctylist[127,1],ctylist[127,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[127,1],ctylist[127,4],ctylist[128,1],ctylist[128,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[128,1],ctylist[128,4],ctylist[129,1],ctylist[129,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[129,1],ctylist[129,4],ctylist[130,1],ctylist[130,4],length=.12,angle=20,col="blue",lwd=3)
arrows(ctylist[130,1],ctylist[130,4],ctylist[131,1],ctylist[131,4],length=.12,angle=20,col="blue",lwd=3)
text(-0.033, 0.08131, "1939", col="red")
text(-0.027, 0.1207, "1940", col="red")
text(0.163, .167, "1941", col="red")
text(0.105, .116, "1942", col="red")
text(0.155, .153, "1943", col="red")
text(0.072, .128, "1944", col="red")
text(0.055,  .0932, "1945", col="red")
title(main="Actor Effects of Russia and the USA, 1816-2001")


#figure 5b
country<-c("USA", "GMY", "UKG", "RUS", "IND", "ITA", "FRN", "CHN", "JPN", "AUS")
ctylist.pe.2 <- data.frame(matrix(NA,length(partner.effect.i),length(country)))				
				for (cty in 1:length(country)){
				for ( i in 1:length(partner.effect.i)){
				ctylist.pe.2[i,cty]<-partner.effect.i[[i]][match(country[cty], row.names(partner.effect.i[[i]]))]		
				}
				}

names(ctylist.pe.2) <- country
ctylist.pe.2[ctylist.pe.2==0]<-NA
years<-1816:2001
	row.names(ctylist.pe.2)<-years

plot(ylim=c(-.1,.20),xlim=c(-.1,.20),bty="n",las=1,ctylist.pe.2[,1],ctylist.pe.2[,4],pch=19,xlab="USA Target Effect",ylab="Russia Target Effect",axes=T,
     col=c(rep("gray67",123),rep("red",length(124:130)),rep("gray50",length(131:174)),rep("gray90",length(175:186))), cex=.9
)
abline(h=summary(ctylist.pe.2[,1])[3])
abline(v=summary(ctylist.pe.2[,4])[3])
arrows(summary(ctylist.pe.2[,1])[3],summary(ctylist.pe.2[,4])[3],ctylist.pe.2[124,1],ctylist.pe.2[124,4],length=.11,angle=20,col="blue",lwd=2)
arrows(ctylist.pe.2[124,1],ctylist.pe.2[124,4],ctylist.pe.2[125,1],ctylist.pe.2[125,4],length=.12,angle=20,col="blue",lwd=2)
arrows(ctylist.pe.2[125,1],ctylist.pe.2[125,4],ctylist.pe.2[126,1],ctylist.pe.2[126,4],length=.12,angle=20,col="blue",lwd=2) #1941
arrows(ctylist.pe.2[126,1],ctylist.pe.2[126,4],ctylist.pe.2[127,1],ctylist.pe.2[127,4],length=.12,angle=20,col="blue",lwd=2) #1942
arrows(ctylist.pe.2[127,1],ctylist.pe.2[127,4],ctylist.pe.2[128,1],ctylist.pe.2[128,4],length=.12,angle=20,col="blue",lwd=2) #43
arrows(ctylist.pe.2[128,1],ctylist.pe.2[128,4],ctylist.pe.2[129,1],ctylist.pe.2[129,4],length=.12,angle=20,col="blue",lwd=2)
arrows(ctylist.pe.2[129,1],ctylist.pe.2[129,4],ctylist.pe.2[130,1],ctylist.pe.2[130,4],length=.12,angle=20,col="blue",lwd=2)
arrows(ctylist.pe.2[130,1],ctylist.pe.2[130,4],ctylist.pe.2[131,1],ctylist.pe.2[131,4],length=.12,angle=20,col="blue",lwd=2)

text(0.036, .128, "1939", col="red", lwd=7)
text(-0.04, .007, "1941", col="red", lwd=7)
text(-.055, -.053, "1942", col="red", lwd=7)
text(-0.016, -0.032, "1943", col="red", lwd=7)
text(-0.0183, -0.049, "1944", col="red", lwd=7)
title(main="Target Effects of Russia and the USA, 1816-2001")

# figure 6
plot(1816:2001, unique.variance, type="l", col="black", main="MIDs and the Social Relation Model", ylab="",xlab="Years",bty="n", axes=FALSE, ann=FALSE,lwd=2)
axis(1, at=c(1815,1914,1939,1960,1991,2001),labels=c("1815","1914","'39","'60","'91","2001"),1,cex=.13)
axis(2, las=1)
text(1968, .04, "Dyadic Variance", col="black",cex=.7)
points(1816:2001, actor.variance, type="l", col="blue")
text(1970, .038, "Actor Variance", col="blue",cex=.7)
points(1816:2001,partner.variance, type="l", col="skyblue")
text(1974, .036, "Partner Variance", col="skyblue",cex=.7)


###########################
### TRADE example (section 4.2)
###########################

#Clean up data a bit
dd<-read.csv(paste(datapath,"data/exIMF2000.CSV",sep="/"), header=T)
rownames(dd)<-dd[,1]
dd<-dd[,2:11]
Y<-as.matrix(dd) 

#run the gbme
gbme(Y,directed=T,k=2,awrite=T,bwrite=T)

out<-read.table("data/OUT",header=T)

#get random effects for senders and receivers 
wg.names<-as.character(rownames(dd))
A<-read.table("data/A",sep=" ",col.names=wg.names)
B<-read.table("data/B",sep=" ",col.names=wg.names)
# throw away the first 5,000 iterations
A<-A[-c(1:50),]
B<-B[-c(1:50),]
a<-apply(A,2,mean) # mean random effect for exporter
b<-apply(B,2,mean) # mean random effect for receiver
cbind(a,b)
cor.test(a,b)

#run the srm on the trade data
outTr<-dyads(Y)

#sort and compare
aiT<-outTr$actor.effect.i[order(outTr$actor.effect.i[,1]),]
aiTgbme<-sort(a)
biT<-outTr$partner.effect.i[order(outTr$partner.effect.i[,1]),]
biTgbme<-sort(b)

# Figure 7 (a)
dotchart(aiT, pch=16, main="SRM Actor Effects for Trade in 2000")
mat2000<-Soc.Matrix$'2000' #from MID data

outW<-dyads(mat2000)
aiW<-(outW$actor.effect.i)

aiW<-outW$actor.effect.i[order(outW$actor.effect.i[,1]),]
aiW<-as.matrix(aiW)
which(rownames(aiW)=="USA") #187
which(rownames(aiW)=="CHN") #184
which(rownames(aiW)=="FRN") # 150
which(rownames(aiW)=="GMY") # 154
which(rownames(aiW)=="ITA") # 160
which(rownames(aiW)=="ROK") # 120
which(rownames(aiW)=="JPN") #121
which(rownames(aiW)=="NTH") # 147
which(rownames(aiW)=="UKG") # 185

namesA<-c("ROK", "JPN", "NTH", "FRN", "GMY", "ITA","CHN", "UKG", "USA")
use<-c(187, 184, 150, 154, 160, 121, 120, 147, 185)
aiW<-aiW[use]
aiW<-sort(aiW)
names(aiW)<-namesA

#figure 7b
dotchart(aiW, pch=16, main="SRM Actor Effects for MID Data in 2000")


# Figure 8 (b)
plot(aiT[1], aiTgbme[1], ylim=c(-6000, 6000), xlim=c(-11000,13000),bty="n",las=1,pch=19,xlab="GBME Effects",ylab="",axes=T)
title(main="Sender Effects for 2000 IMF Trade Data")
title(ylab="SRM Effects", line=4.2)
points(aiT[1], aiTgbme[1], pch=19)
text(-9300, -4900, "S. Korea")
points(aiT[2], aiTgbme[2], pch=19)
text(-7600, -2600, "Italy")
points(aiT[3], aiTgbme[4], pch=19)
text(-6668, -1280, "Netherlands")
points(aiT[4], aiTgbme[3], pch=19)
text(-2500, -3050, "Hong Kong", pch=19)
points(aiT[5], aiTgbme[5], pch=19)
text(-4204, -500, "UK")
points(aiT[6], aiTgbme[7], pch=19)
text(-2700, 370, "France")
points(aiT[7], aiTgbme[6], pch=19)
text(-90, -1400, "China")
points(aiT[8], aiTgbme[10], pch=19)
text(10371, 6050, "USA")
points(aiT[9], aiTgbme[9], pch=19)
text(11500, 4523, "Germany")
points(aiT[10], aiTgbme[8], pch=19)
text(12400, 200, "Japan")

# Figure 8 (b)
par(mar=c(5,5,4,2) + 0.1)
plot(biT[1], biTgbme[1], xlim=c(-15000, 30500), ylim=c(-11000,15000),bty="n",las=1,pch=19,xlab="GBME Effects", ylab="",axes=T)
title(main="Target Effects for 2000 IMF Trade Data")
title(ylab="SRM Effects", line=4.2)
points(biT[1], biTgbme[1], pch=19)
text(-11100, -10000, "S.Korea")
points(biT[2], biTgbme[2], pch=19)
text(-1900, -7300, "Hong Kong")
points(biT[3], biTgbme[4], pch=19)
text(-7700, -2700, "Italy")
points(biT[4], biTgbme[3], pch=19)
text(-100, -5500, "Netherlands")
points(biT[5], biTgbme[6], pch=19)
text(-4660, 200, "China")
points(biT[6], biTgbme[5], pch=19)
text(3500, -3775, "Japan")
points(biT[7], biTgbme[8], pch=19)
text(-1070, 2401, "France")
points(biT[8], biTgbme[7], pch=19)
text(1900, -100, "UK")
points(biT[9], biTgbme[9], pch=19)
text(7808, 7144, "Germany")
points(biT[10], biTgbme[10], pch=19)
text(30000, 14700.71, "USA")


###########################
### MORROW EXAMPLE (4.3)
###########################

# SRM analysis
Y<-morrow.pow5.matrix
out<-dyads(Y)

# basic ranked plot
ncty<- dim(Y)[[1]]
sort.cty<-out$actor.effect.i[order(out$actor.effect.i[,1]),]
sort.cty.pa<-out$partner.effect.i[order(out$partner.effect.i[,1]),]

cor(out$actor.effect.i[,1], out$partner.effect.i[,1])
cty.name<-morrow.pows.t5.dem$violatorfullname

#figure 9 (a and b)
names.sorted.ai<-c("Cyprus", "Turkey", "China", "Vietnam", "Argentina", "Saudi Arabia & allies", "United Kingdom", "Tanzania", "Syria", "Israel", "Uganda & Libya", "USA & allies", "Iran", "Ethiopia & Cuba", "Somalia", "Iraq")
names.sorted.pa<-c("Cyprus", "Turkey", "China", "Vietnam", "United Kingdom", "Argentina", "Saudi Arabia & allies", "Tanzania", "Syria", "Israel", "Uganda & Libya", "USA & allies", "Iran", "Somalia", "Ethiopia & Cuba",  "Iraq")
dotchart(sort.cty, labels=c(names.sorted.ai), pch=19) 
dotchart(sort.cty.pa, labels=c(names.sorted.pa), pch=19)

#figure 10
sort.cty.uniq<-out$actor.effect.i[order(out$unique.effect.ij[,1]),]
sort.cty.uniq.bi<-out.bi$actor.effect.i[order(out.bi$unique.effect.ij[,1]),]
uniq.net<-as.matrix(out$unique.effect.ij)
diag(uniq.net)<-0
plot.network(network(uniq.net, directed=T), label=cty.name, usearrows=T,edge.col=8, vertex.col="darkblue",
label.col="black", label.pos=1, label.cex=.75, edge.lwd=.1)
library(qgraph)
cty.name<-morrow.pows.t5.dem$violatorfullname
cty.name.abb<-list("ARG", "UK", "CYP", "TAZ", "SOM", "IRN", " TUR", " IRQ", "SYR" , "ISR", "CHN", "DRV", "ETH", "UGA", "USA", "SAU")
qgraph(uniq.net, minimum=-0.7619048, maximum=  3.6041667, cut=0, labels=cty.name.abb, asize= .1, arrows=FALSE, label.scale= TRUE)
par(fig=c(.8,1,.8,1),new=T)
qgraph(morrow.pow5.matrix, weighted=TRUE, labels=cty.name.abb)
par(fig=c(0,1,0,1))


#GBME (code to run the GBME)
Y<-morrow.pow5.matrix
n<-dim(morrow.pow5.matrix)[1]
Xs<-as.matrix(morrow.pows.t5.dem$violator_democracy7)
X.array<-array(dim=c(n,n,2))
X.array[,,1]<-morrow.war5.matrix
X.array[,,2]<-morrow.ally5.matrix
Xd<-X.array

gbme(Y=Y,Xd=Xd,Xs=Xs,fam="gaussian",k=2, awrite=T, bwrite=T, ofilename="morrow.pow5.out", 
efilename="morrow.pow5.U", ffilename="morrow.pow5.V", afilename="morrow.pow5.A", 
bfilename="morrow.pow5.B", NS=500000, odens=100)

# Loading in GBME results
morrow.pow5.out<-read.table(file="data/morrow.pow5.out", header=T, sep=" ")
st<-3001
fin<-5000
morrow.pow5.U<-read.table("data/morrow.pow5.U", header=F)
morrow.pow5.V<-read.table("data/morrow.pow5.V", header=F)

# table 3 (replication of MORROW's ordered probit)
load("data/newdata3.r") #dataframe for OLS, o probit, and logit
model.oprobit <-polr(as.factor(Vio_comply) ~ Vic_comply + dem7_sender + ally + morrow.pow5.b, data=newdata3, Hess=TRUE, model=TRUE,method=c("probit"))
#table
summary(model.oprobit)

#Table 4 (estimates and credible intervals )
constant<-morrow.pow5.out$b0[3001:5000]
mean(morrow.pow5.out$b0[3001:5000]) #constant
quantile(constant,probs=0.025)
quantile(constant,probs=0.975)

mean(morrow.pow5.out$bs1[3001:5000]) # democracy-sender level
quantile(morrow.pow5.out$bs1[3001:5000],probs=0.025)
quantile(morrow.pow5.out$bs1[3001:5000],probs=0.975)

mean(morrow.pow5.out$rho[3001:5000])  # reciprocity
quantile(morrow.pow5.out$rho[3001:5000],probs=0.025)
quantile(morrow.pow5.out$rho[3001:5000],probs=0.975)

mean(morrow.pow5.out$bd1[3001:5000]) # at war- dyadic level
quantile(morrow.pow5.out$bd1[3001:5000],probs=0.025)
quantile(morrow.pow5.out$bd1[3001:5000],probs=0.975)

mean(morrow.pow5.out$bd2[3001:5000]) # alliances- dyadic level
quantile(morrow.pow5.out$bd2[3001:5000],probs=0.025)
quantile(morrow.pow5.out$bd2[3001:5000],probs=0.975)

mean(morrow.pow5.out$s2a[3001:5000]) # actor variance/common sender
quantile(morrow.pow5.out$s2a[3001:5000],probs=0.025)
quantile(morrow.pow5.out$s2a[3001:5000],probs=0.975)

mean(morrow.pow5.out$sab[3001:5000]) # sender-receiver
quantile(morrow.pow5.out$sab[3001:5000],probs=0.025)
quantile(morrow.pow5.out$sab[3001:5000],probs=0.975)

mean(morrow.pow5.out$s2b[3001:5000]) # receiver
quantile(morrow.pow5.out$s2b[3001:5000],probs=0.025)
quantile(morrow.pow5.out$s2b[3001:5000],probs=0.975)

mean(morrow.pow5.out$s2e)



